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Abstract 



ABSTRACT We propose a dissipative electro-elastic network model (DENM) to describe the 
dynamics and statistics of electrostatic fluctuations at active sites of proteins. The model com- 
bines the harmonic network of residue beads with overdamped dynamics of the normal modes 
of the network characterized by two friction coefficients. The electrostatic component is intro- 
duced to the model through atomic charges of the protein force field. The overall effect of the 
electrostatic fluctuations of the network is recorded through the frequency-dependent response 
functions of the electrostatic potential and electric field at the active site. We also consider the 
dynamics of displacements of individual residues in the network and the dynamics of distances 
between pairs of residues. The model is tested against loss spectra of residue displacements and 
the electrostatic potential and electric field at the heme's iron from all-atom molecular dynamics 
simulations of three hydrated globular proteins. 

Key words: protein electrostatics; dissipative dynamics; elastic network; electrostatic re- 
sponse function; loss spectrum 
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INTRODUCTION 

Proteins in solution exist not as static structures but instead as dynamic entities interconverting 
between conformational sub- states on the time- scale reflecting the activation barriers involved 
in the transitions It is becoming increasingly clear that biomolecules are fluctuating ma- 

chines, using dynamical equilibria between their conformational states to promote their function 

m 

. The flexibility of the protein structure naturally leads to changes in both the charge distri- 
bution of the protein itself and in the polarization of the interfacial waters following the protein 
motions. These dynamically fluctuating polarization fields will affect many properties sensitive 
to electrostatics, including the rates of enzymatic reactions (9). The dynamic nature of proteins, 
enhanced by the conformational manifold provided by hydration water (1), requires the shift 
of the focus from statistical averages and corresponding thermodynamic parameters, such as 
equilibrium Gibbs energies, to entire statistical distributions of observables and their dynamics 
expressed through time correlation functions. 

This requirement is certainly the case for the problem of protein electron transfer, where, 
according to the standard Marcus picture (@), both the average of the energy gap between the 
donor and acceptor energy levels of the electron and its fluctuation are required to determine 
the probability of electron tunneling. This is just one example when the knowledge of the entire 
fluctuation spectrum of the observables is of significant interest [l). Other applications, such 
as electrostatic fluctuations reported by broad-band dielectric spectroscopy ( 8), THz absorption 
d^, and light scattering techniques ([Iq) require either the entire fluctuation spectrum or a set of 
cumulants of the corresponding experimental observables. 

With the general focus on the fluctuation spectrum of proteins, we propose here a model for 
the calculation of the statistics and dynamics of the electrostatic potential and electric field in 
proteins. Electrostatic interactions are clearly important for protein stability and function (J). 
Electrostatic solvation and interactions affect the stability of folded proteins and their aggrega- 



tion and crystallization (Illl4l3). Electrostatics might also play a significant role in promoting 
high-temperature flexibility of proteins through the coupling of the charge distribution in the 
protein to the electrostatic fluctuations produced by the protein- water interface 

The dynamics of electrostatic fluctuations spans an enormous range of time-scales from 
sub-picoseconds to sub-microseconds, and possibly longer ([si 15, 16i). While spectroscopic 
techniques are capable of recording the ultra-fast fs-to-ps relaxation (Il7|) . the fluctuations of the 
protein-water interface recorded by Mossbauer spectroscopy and neutron scattering cover much 
longer time-scales from sub-nanoseconds to sub-microseconds (3). Our present model aims at 
these long, and perhaps even longer, time-scales by coarse-graining the protein into an elastic 
network of beads, each representing a protein residue. 

Elastic network models have consistently shown good performance in characterizing global. 



large-scale motions of proteins (Il8l- l27|) . While the low-frequency portion of the spectrum of 



protein motions is mostly determined by the distribution of mass and molecular shape (l28l-l30l). 



improvements are still needed to account for deficiencies of networks when localized protein 



motions are involved (13 11-1341) . Electrostatic interactions, on the other hand, are notoriously 
long-ranged. The electrostatic potential, slowly changing over a nanometer- size biomolecule, 
effectively averages out local structural variations and is mostly sensitive to the global distri- 
bution of charge within the molecule. From this perspective, the combination of even a basic 
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elastic network with the molecular charge distribution might be sufficient to describe the statis- 
tics of the potential fluctuations and their slow dynamics. 

We propose here a model combining the distribution of molecular charge from the standard 
atomic force-fields with the dynamics and statistics of protein motions derived from an elas- 
tic network. To test the model, we employ all-atom Molecular Dynamics (MD) simulations 
of three hydrated heme proteins. Previous analysis of these data has emphasized the impor- 
tance of nanosecond (ns) relaxation modes in the electrostatic fluctuations of the protein ma- 
trix and the protein-water interface (Il4ll35|) . This time-scale is relevant to biological function 
since heme proteins are typical components of biology's energy chains transporting electrons 
on nanosecond-to-microsecond time-scales (^Jsd). The nanosecond time-scale fluctuations are 
also consistent with the dynamics of the global motions of the network of residues. An elastic 
network model naturally fits the physics of the problem. The present contribution is a first step 
in developing network models of the electrostatic response of hydrated proteins. Here, we fo- 
cus on the electrostatics of the protein matrix only, leaving the water component of the overall 
response to future studies. 



THEORY AND METHODS 



Elastic network model 

A coarse-grained elastic network model (ENM) assigns a node to a collection of atoms reducing the 
computational burden of an all-atom normal-mode analysis. The typical coarse-graining is done on the 
level of individual aminoacids (1-bead model (l37h). The position of the node is defined by the coordinates 



of the C" atom. The springs connecting the nodes represent the bonded and non-bonded interactions 



within an accepted cutoff distance (|18M20t) or by a distance-dependent force constant 
The ENM diagonalizes the Hessian matrix H derived from a simplified Hookean potential suggested by 
Tirion (flsh . This potential, Ejj = C{rij — ro,,y)^/2, describes the elongation rjj = |r,- — rj] between nodes 
/ and j in the network characterized by one universal force constant C and the structural information 
stored in the equilibrium bead positions ro.,- {rojj = \rQj — rQj\). 

The Hessian Hjj is a 3N x 3N matrix representing the protein elastic energy as a quadratic form of 
the Cartesian displacements 5rf = rf — r^- 



of individual beads in the network 



E = {C/2)Y,Hf5rf5/j 



(1) 



Here, / and j run between 1 and A^, a, P indicate the Cartesian projections, and summation over repeated 
Greek indices is assumed here and throughout. The Hessian is diagonalized by the unitary matrix U 
producing 3N eigenvalues km- This standard linear algebra formalism (jlSL Il9h yields the statistical 
correlator between the displacements of residues / and j in the network 



(8rf 5r5) = (PC) 



-1 



m.y 



(2) 



where P = 1/ {k^T) is the inverse temperature. 

If the network is characterized by a universal force constant and a cutoff, it yields a bell-shaped 
density of vibrational states. A refinement of that network by adopting a stronger coupling between 
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covalently bound neighbors splits the density of states into two maxima, in better agreement with all- 
atom calculations (24, 4^ 41). This approximation is adopted in our present calculations: the spring 
constant was multiplied by a constant factor £ = 100 for covalent neighbors. In addition, a uniform 
cut-off radius of 15 A was used in all calculations. 



Overdamped network dynamics 



The standard mechanical ENM outlined above obviously lacks dissipative dynamics. The equations 
of motion are harmonic, implying oscillatory time correlation functions. In contrast, most correlation 
functions of hydrated proteins observed by scattering |42|) and relaxation (|43|) techniques are 

exponential, corresponding to the overdamped (Debye) dynamics, or stretched-exponential (8). Several 
approaches can be implemented to incorporate dissipative dynamics into the mechanical network of 
beads. Langevin equations of motion for the E NM p otential, within the general framework of the Lamm- 
Szabo formalism (44), have been suggested (|34, 48). This approach, and some early suggestions (45, 46), 
still requires parameterization (done by fitting the rotational and translational diffusion coefficients from 
MD) and, in addition, doubles the size of the matrix to be diagonalized. We have adopted here a more 
straightforward formalism not requiring additional computational resources. 

Instead of a harmonically oscillating network, we have assumed for each normal mode q„„ diago- 
nalizing the network's Hessian, an overdamped motion described by the dissipative memory kernel C,{t) 
(l47h . The equation of motion for such overdamped dynamics is 



f i;{t - t')q,n{t')dt' + 'k,nqm = 'P{t), 
Jo 



(3) 



where F{t) is an external force. 

The distinction between Eq. ^ and the Langevin network (|34il48|) is worth emphasizing here. In the 
latter, uncoupled Langevin dissipative equations are first assigned to each bead of the network. However, 
it was noted that dissipative equations for the beads are likely to become coupled when the Langevin 
dynamics of beads are consistently derived by integrating out the fast degrees of freedom in the equations 
of motion (49). Given that the assignment of uncoupled Langevin equations to the individual beads is 
phenomenological from the onset, one can introduce similar phenomenology at a different level of the 
theory. Here, in the spirit of the standard normal-mode analysis, we introduce uncoupled dissipative 
equations for the normal modes (Eq. (|3]l). This is still a phenomenological assumption requiring further 
testing, but it reduces to the standard normal-mode analysis in the static limit. 

When Laplace-Fourier transform (|47h is applied to Eq. the displacement response function for 
collective mode q,„ follows 



Xm(CO) 



(4) 



where ^((o) is the Laplace-Fourier transform of the friction kernel ^{t). Extending this procedure to all 
normal modes of the network, one obtains the response function of bead displacements 



(5) 



In this equation, the response function x,y((o), which is a rank-2 tensor, represents the displacement of 
residue / of the protein due to an oscillatory force with frequency co applied to residue j and propagated 
through the network to residue /. This response function is the basis of the dissipative electro-elastic 
network model (DENM) proposed here. At co = 0, Eq. dD returns the standard result of the ENM given 
by Eq. dUl. 
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active site 



Figure 1 : Cartoon illustrating the application of the dissipative electro-elastic network model (DENM) 
to the calculation of the electro-elastic response function at an atomic position in the active site of the 
protein. The elastic displacements of residues / and j create fluctuations of the electric fields by atomic 
charges qi^ and qj^. Electric fields E,jt and Ej^ are produced by the con^esponding charges qnc and qj^ at 
the active site. The overall electric field of residue / at the active site, Eq,-, is obtained by summation of 
all E,jt produced by the charges qn^ of the residue. 



Electrostatic potential response function 

The ENM has shown good performance for low-frequency structural/mechanical deformations of indi- 
vidual proteins and biomolecular assemblies (l27h . Our goal here is to supplement these low-frequency 
motions with atomic partial charges to model the dynamics and statistics of electrostatic fluctuations 
at active sites of proteins. Coarse-grained electrostatics of proteins has been addressed in a number of 
recent papers (50-52), mostly dealing with long-range protein assembly and interactions (|53). These ap- 
proaches also involve coarse-graining of the charge distribution (15 ll. 1521) or generating effective charges 
to reproduce the protein's external field Given that we are interested in the electrostatic fluctua- 
tions, which are more sensitive to a local charge distribution, atomic charges from force-field potentials 
are more consistent with our purpose. 

Our approach starts with solving the standard ENM problem to obtain a set of displacements 6r,>„ 
for each normal mode m with the eigenvalue X„,. Each of these displacements is then applied to all 
charges qik of residue /, where index k runs over all atoms of residue / (Fig. [T]). Displacements 5r,„, in 
turn produce dipole moments 5/i|^'' = ^i/tSr™ which interact with charges of the protein active site. The 
overall electrostatic effect of the fluctuating protein medium is given as a sum over all contributions from 
each partial charge of all the residues except for the active site. This cumulative effect of all protein 
charges is given by the frequency-dependent response function of the electrostatic potential (co) 



X^((o) = -^£o"yXy(co)4- (6) 



Here, Eq^ is the a Cartesian projection of the electrostatic field produced by all charges qik of residue / at 
an atom in the active site (Fig. [T]) 

Here, tq is the position of the atom in the active site and rtk is the position of charge qik of residue /. 
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Electrostatic field response function 

The potential response function X(^{(i)) in Eq. ^ describes an oscillating electrostatic potential pro- 
duced by the protein matrix in response to placing an oscillatory charge at the position in the active site 
where the potential is recorded (Fe ion in our MD simulations). Similarly, the response function of the 
electric field considers the field produced by the protein in response to an oscillating dipole moment 
mo(0 = mo (co) exp [/cof] placed in the active site. The deformation of the protein matrix induced by this 
probe dipole results in the electric field Eo((o) at the same site. It can be found by summation over the 
contributions from all individual dipoles 5/i,/t((o) arising from residues' atomic charges 

Eo{(o) = Y,Tif5niti(o). (8) 

ik 

Here, T,x. is the dipolar tensor connecting the position of the active site with the charge qn; of residue /. 
The dipole moment 5/x,i:((o) at the position of qik is caused by the residue displacement 5r,i(-((o) and can 
be expressed in terms of the displacement response function 

5fj,ik (ro) = <?i,t £ Xij ■ 8F^- (to) , (9) 

j 

where 6Fy((o) = T,kQjkT^jk ■ nio((o) is the force caused by the dipole mo((o) at residue j. We therefore 
obtain a linear relation between the dipole m()((o) placed at the active site and the electric field Eo((o) 
produced by the protein in response to this perturbation. The proportionality coefficient between the 
external perturbation and the response is the response function given by the rank-2 tensor 

xf (CO) = I ,,.7;fxJ(a))r,?^;7. (10) 

i:j,kJ 

Here, as above, the summation is done over the repeated Greek indices denoting the Cartesian compo- 
nents of the corresponding tensors. We will be mostly interested in the trace of the tensor 

X£N=%rN- (11) 
Response and correlation functions 

The response functions calculated by the DENM formalism need to be related to time correlation func- 
tions supplied by the MD trajectories. The dynamics of a general dynamical variable bX (t) = X {t) — {X) 
is characterized by the time self-correlation function 

Cx{t) = {dX{t)bX{0)). (12) 

The normalized correlation function Sx{t) = Cx{t)/ {{bX)'^) is related to the response function by the 
fluctuation-dissipation equation which forms the basis of our analysis (55,) 

Xxiz) = mxf)[^ + izSx{z)], (13) 

where Sx{z) is the Laplace-Fourier transform of S{t) defined in the upper half of the complex plane of z. 

Since our main focus is on variances of physical properties, we will define the generalized compli- 
ance for the variable X as follows 



/»oo 

Xx = P((5X)2)/2= / Xx(co)(^a)/7ra)). 
Jo 



(14) 
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The same property can be calculated from the co = value of the real part of the response function 

'kx = {\/2)y!x{0). (15) 

In case of residue displacements considered as variable X {X = r), X,- <xC^ is proportional to the inverse 
force constant of the network springs. As for the compliance of a macroscopic body, defined as the 
inverse of stiffness, X,. depends on the shape and boundary conditions, in contrast to elastic moduli 
representing material properties. 

For the elastic network, the generalized compliance reports on how the softness of the network affects 
the variance of the variable of interest. In the case of X = (j) reporting on the potential ([) at the position of 
the heme's iron, X = e^X^ is the reorganization energy of a half redox reaction corresponding to changing 
the oxidation state of the heme (6, 56). The standard definition of half -reaction reorganization energy 
involves, instead of one atom, the distribution of the electronic density of the transferred electron over a 
few atoms of the active site. We simplify this problem here by assuming all electron charge e localized 
at the centroid of this charge density, the iron atom of the heme. The electrostatic potential at a single 
atom is a well defined physical property, but its variance is only an approximate representation of the 
observable reorganization energy of changing the redox state Q). 

We want to gain insight into the dissipative dynamics of the elastic network and its comparison with 
the dynamics calculated from MD simulations. The loss function Xx(®) provides direct access to both 
the set of characteristic relaxation frequencies and their relative contributions. However, it does not 
weigh the low and high frequencies as they contribute to the variance in Eq. (fT4)) . Therefore, in addition 
to the loss function, we will consider the function 

71(0 x^(0) 

This function is normalized, ax{(i))d(i) = 1, and tends to the characteristic relaxation time (x) in the 
CO hmit 

limax((o) = (2/7i)(x), (17) 

where 

poo 

(x) = / Sx{t)dt. (18) 







Proteins used in case studies and the simulation protocol 

Three heme proteins, oxidized form of myoglobin (metmyoglobin, metMB, lYMB) reduced state of 
cytochrome B562 (cytB, 256B), and reduced state of bovine heart cytochrome c (cytC, 2B4Z) were 
simulated by all-atom MD. The parameters of the proteins were taken from CHARMM27 force field 
(57), and NAMD (|58j) was used for the trajectories production. Each protein was placed at the center 



of a cubic box with the side length of ~ 108 A and solvated with TIP3P waters (15 9h . The number of 
waters used in simulations were: 32891 (metMB), 33268 (cytB), and 33189 (cytC). The VMD (6d) 
plugin Autoionize was used to neutralize the simulation cell by adding Na+ and Cl^ ions to bring the 
ionic strength of the solution to 0.1. Particle-mesh Ewald with the grid resolution < 1 A was used 
for electrostatic interactions, and all other non-bonded interactions were calculated within 12 A cutoff. 
Following energy minimization, each protein was simulated for 5 ns in a NPT ensemble atP = I atm and 
T = 300 K. Temperature and pressure were controlled by using the Langevin dynamics with the damping 
coefficient of 5 ps^^. The production NVE trajectory was started at the end of each 5 ns trajectory. 
The NVE ensemble is required for the proper sampling of the long-time dynamics since we found that 
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residue No. residue No. residue No. 

Figure 2: Residue mean-square displacements (msd's) ((5r,)^) (a) and variance-covariance matrix (5r,- • 
bvj) of residues' C" carbons calculated from MD (b) and from DENM (c). The diagonals in (b) and (c) 
represents the residue msd's ((6r,)^). The network force constant in the DENM calculations is k^T /C = 
1 A^. The DENM values indicated by red squares in (a) need to be multiplied by about a factor of two in 
order to obtain the best-fit agreement with the MD data. They are separated for a better visibility in the 
plot. 



using NPT and NVT ensembles artificially accelerates the observables relaxing on sub-nanosecond to 
nanosecond time-scale ( |6l|) . The integration step was 2 fs and simulation frames for the analysis were 
saved each 0.05 ps. The simulation trajectories were 65 ns (metMB), 100 ns (cytC), and 123 ns (cytB). 

The force field standard parameters of the heme in the reduced state are taken from CHARMM27 
(1570 . The charge of iron is ^pe = 0.24 in the reduced state. No standard parameters are available for the 
oxidized heme required for the simulations of metmyoglobin. These force field parameters were taken 
from Autenrieth et al ( |62|) . The atomic charge of oxidized Fe is = 1.34 in this parametrization. The 
iron atom is five-coordinate in its oxidized state. In order to allow the iron to shift out of the porphyrin 
plane, the harmonic force constant of the bond between Fe and Ng2 of histidine 93 was adjusted to 65 
kcal/mol as suggested in Ref. (,631) and also implemented in Ref. (iM)- 



RESULTS AND DISCUSSION 

Statistics and dynamics of residue displacements 

Consistent with many previous studies, we have found that elastic networks are capable of re- 
producing the basic pattern of the distribution of mean-square displacements (msd's) along the 
protein backbone. The left panel in Fig.[2lcompares residue msd's from MD with DENM calcu- 
lations, while two other panels show the maps of variance-covariance matrices. Overall, there 
is a good agreement between the alteration in residue displacements along the backbone, and 
corresponding cross-correlations, calculated from the two sets of data. However, the network 
force constant required for the best fit of the msd from MD simulations, c:^ 0.3 kcal/(mol A^), 
is somewhat lower than the value of ~ 1 kcal/(mol A^) typically adopted from fitting the B- 



factors from crystallography (| 181-120). This value is also about a factor of two lower than 0.6 
kcal/(mol A^) adopted below to globally fit the variances of the field and electrostatic potential 
fluctuations from MD data for all three proteins studied here. 

The fitting of the force constant to the absolute msd's from MD makes the network too soft. 
This outcome is expected since the elastic network cuts high-frequency vibrations from the 
density of states and thus underestimates the absolute magnitudes of the residue displacements 



(|27h . It appears more consistent to use for the purpose of fitting only the portion of the msd 
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co/ns co/ns 

Figure 3: Loss spectra x'l {(iij/x'iiO) of individual residues (/ = 21, a) and of pairs of residues 
X'i'j{0)) /x'ij{0) (/ = 21 and j = 84, b). The results have been obtained from MD (solid hnes) and DENM 
(dashed lines). The calculations are done for cytB with the residues used in the calculations shown in Fig. 
m The DENM calculations are done with Q = 30 ns used for the electrostatic calculations in Fig. [5] (red 
dashed line) and with ^/ = 5 ns (blue dash-dotted line). The rest of DENM parameters are: = 0.006Q, 
a = 0.35, 8 = 100. 



related to protein's global motions. This component can in fact be separated from the vibrational 
part in the temperature dependence of the msd since global motions of the protein become 



observable at high temperatures, above the temperature of the protein dynamical transition (65). 
This high-temperature portion of the protein msd is about a half of the total magnitude (66% 
which is close to the factor required to reconcile the force constants from electrostatic variances 
and msd's. 

Figure[3]shows the dynamics of positions r, (?) of individual residues within the network and 
the dynamics of distances between residues rij{t) = r,(?) —rj(t). It reports the loss functions 
Xf((o) and %|y((0) calculated from the self-correlation functions of individual residues 

Q(0 = (5r,(0-6r,(0)) (19) 

and self-correlation functions of distances between the residues 



Qjit) = {dr,j{t)-r,jiO)). (20) 

These latter correlation functions are experimentally available from FRET (67) recording the 
evolution of the distance between two residues tagged with chromophores. 

The comparison of xf((0) from DENM and MD is illustrated in Fig. [3^ for only one residue, 
i — 2l, belonging to cytB protein. The pair of residues with / = 21 and j — 84 is used to calculate 
Cij{t) in Fig. [3h. The position of the pair in the backbone of cytB is shown in Fig. ID Examples 
of calculations for several additional pairs of residues from cytB can be found in Supporting 
Materials. 

Both correlation functions, C,(?) and C,j(?), report the dynamics on roughly two time-scales, 
seen as two peaks in the their loss functions in Fig. [3l The two peaks represent, with varying 
weights, the faster local dynamics of individual residues and the slower global dynamics of the 
network. 

The elastic network with single-exponential, Debye dynamics does not capture two time- 
scales of the dynamics from MD simulations. When the Debye memory kernel ^((o) = 
is used in the response function %m(to) in Eqs. (H)) and ([5]), the loss function shows only one 
relaxation peak, slightly deformed from a simple Debye form by the distribution of network's 
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Figure 4: Cartoon of cytB showing tlie residues (red) used to produce tlie correlation functions C,(f) 
(/ = 21) and Cij{t) = 21,84, Eqs. ( [T9l ) and (l20l) ) and corresponding loss spectra X-'((o) and X;y(w) 
shown in Fig. [3l The residues used in producing xJy((o) are connected by the red soUd line. The red 
spheres representing the two chosen residues are centered at their corresponding C" carbons. The Fe 
atom of the heme is rendered as a yellow sphere. 



eigenvalues X,,,. A single friction coefficient evidently misses the fact that energy dissipation 
decreases with increasing frequency of vibrations. Since the protein network is characterized by 
two spring constants, lower for non-covalent neighbors within the cutoff and higher for covalent 
neighbors, it must possess at least two characteristic friction coefficients. 

Friction kernels typically used in applications are phenomenological /^v") and we as well 
cannot offer a consistent theoretical formalism capturing the complex dynamics of residue dis- 
placements. We have therefore adopted a phenomenological response function best describing 
the MD results. Instead of searching for a functional form of ^(co) reproducing MD simulations, 
we have resorted to a two-Debye form for the entire response function of bead displacements in 
Eq. © 

/ \ « I— a ^^^^ 

Xm(0)) = . ^ , . + . r , 1 • (21) 

Here, the amplitude a represents the relative weight of the high-frequency dissipation with the 
high-frequency friction ^/,. Correspondingly, ^/ is the low-frequency friction, which is larger in 
magnitude. It is obvious that Eq. (|2T)) satisfies the static limit %m(0) =Xyy^ . 

The two-Debye form of is directly applied to calculate the response function of 
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residue displacements 

(22) 

where the function x"'^(a)) is given by Eq. Q. The loss spectra of residue displacements show 
two peaks and qualitatively agree with the simulations. The positions and relative heights of the 
peaks can be adjusted by choosing C,h i and the amplitude a in Eq. (|2TI) . as is illustrated in Fig. [3^ 
by dashed and dash-dotted lines. The model can therefore be well parameterized to reproduce 
the dynamics of displacements of a small subset of residues. However, it fails to discriminate 
between the dynamics of residues across the protein backbone. 

The loss functions calculated from MD for a number of residues (see Supporting Materials) 
display similar two-peaks pattern, but the weights and positions of their peaks vary among 
the residues. In contrast, the dynamics of displacements of individual residues in the elastic 
network are mostly driven by the global motions of the entire network. As a result, there is 
little difference between the loss functions of individual residues calculated with DENM. The 
same statement applies to the dynamics of distances between residues shown in Fig.[3h. While 
one can reproduce the loss function of a chosen pair of residues by adjusting the parameters of 
Xm(ti>) in Eq. (|2T]) . these parameters do not translate to all pairs in the network and instead can 
be only viewed as an average representation of the pairs dynamics. 

The fast component of the dynamics, prominent in the loss function %^ ((o) , is less significant 
in the function ax(tt)) due to the l/co scaling of X^((o) in Eq. (fT6l) . The function ax((o) is 
more relevant to problems related to the generalized compliance as defined by Eq. (fT4l) . The 
correct representation of the fast dynamics is therefore of lesser importance for these type of 
problems. The l/co scaling difference between x'xi^) ^^'^ 0(x(«>) explains the relative success 
of the network models in reproducing the pattern of residue msd's (Fig.O. Those are obtained 
by frequency integration of a^((o) and are less sensitive to the details of local dynamics of 
individual beads in the network represented by high-frequency peak of their loss functions. 

Dynamics of the electrostatic field and potential 

Essentially all time correlation functions obtained here from MD simulations show a pattern that 
is roughly represented by a two-exponential decay, requiring in some cases a third decaying ex- 
ponent for a better mathematical fit. The relative weights of the fast and slow components of the 
correlation functions differ depending on the observable. The loss function of the electric field 
%^((jo) at the position of protein's heme is mostly single-exponential, with the low-frequency 
peak much exceeding the high-frequency one. This pattern is less uniform for the loss function 
of the electrostatic potential showing different weights of the low- and high-frequency peaks 
in %^'((o) among the proteins dssj). Figure [5] shows the loss spectra x^((o)/x^(0) calculated 
from the DENM and MD for cytB and metMB proteins. Similarly to the case with the residue 
displacements, the low-frequency coefficient needs adjustment to reproduce the position of the 
main loss peak: the value of Q = 30 ns was adopted for cytB and C,j = 10 ns was taken for 
metMB. Still, the set of parameters taken to reproduce the loss spectrum of electrostatic poten- 
tial does not perform as well when applied to the electric field (cf. Figs. [5]: and[5ll). 

The relative contribution of the fast dynamics component is reduced in the spectral function 
(Xxioo) describing the frequency variation of the generalized compliance (Eqs. (fT4l and (fT6l)). 
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Figure 5: Loss spectra x'i^{oi)/x'^{0) (a,c) and j!e{'^)/j!e{^) (b,d) for cytB (a,b) and metMB (c,d). The 
results are from MD trajectories (dashed lines) and from the DENM calculations (solid lines). The 
DENM calculations were done with the two-Debye Xm(ro) in Eq. (|2TI ). The two-Debye relaxation pa- 
rameters are: ^, = 30 ns, = 0.006Q, a = 0.35 for cytB and ^, = 10 ns, = 0.0002^,, a = 0.35 for 
metMB. The elastic network is defined with k^T /C = \ K^,z = 100, and the cutoff radius of 15 A. 
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Figure 6: Functions ax(G)) (ns) for cytB, cytC, and metMB from MD (dashed lines) and DENM (sohd 
lines). The results for the electrostatic potential (a,|,((o), blue) and electric filed (a£(a)), red) are shown. 
A single set of network parameters is used for all three heme proteins: ^/ = 30 ns, = 0.006Q, a = 0.35, 
8 = 100. 



Figure [6] compares ax(co) from DENM with MD data for all three globular proteins studied 
here. In order to show the ability of DENM to describe the entire set of MD data, a single 
set of network and friction parameters has been assigned to all three proteins. As expected, 
the performance of DENM is better for ax((0). The average relaxation time (x), given by the 
CO —J- limit of CXx((o) (Eq. (flTl)). is well reproduced by the network calculations. The deviations 
for some observables correspond to the cases when the dynamics are strongly dominated by its 
high-frequency component. For these cases, the fit can be improved by adjusting the weight a 
of the fast component in Eq. (|2TI) to a higher value. 

The loss functions Xx{^) shown in Figs. [3l [51 and [6] are normalized with the corresponding 
values of X^(0) and therefore are not affected by the magnitude of the spring constant C de- 
scribing the non-covalent interactions within the cutoff distance = 15 A. The absolute values 
of the variances are, however, proportional to k^T /C, as for instance in Eq. We adopted in 
our calculations the value of k^T /C = 1 A^, which is equivalent to C ~ 0.6 kcal/(mol A^). The 
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Table 1 : Generalized compliances (Eq. (fT4l) ) for the electrostatic potential fluctuations Xi}^ (also 
known as reorganization energies of redox half reaction X = e^X,^, kcal/mol) and for the electric 
force, F = eE (kf, kcal/(mol A^)) for the three proteins used in the case study; the network is 
assigned the force constant of C = 0.6 kcal/(mol A^) for non-covalent springs and £C, £ = 100 
for covalently bound neighbors. The cutoff radius is 15 A. 



Protein 


e% (DENM) 


e% (MD) 


If (DENM) 


If (MD) 


MetMB 


85 


228 


436 


243 


CytB 


115 


136 


9 


17 


CytC 


141 


159 


18 


83 



results of calculating the generalized compliances Xx for electrostatic observables (Eqs. (fT4l) 
and (fTSi) ) with this choice of the force constant are summarized in Table [U The agreement is 
semi-quantitative at best. While the value of the force constant can clearly be adjusted in each 
particular case, this change does not propagate into equally good agreement between X,^ and Xf, 
where F = eE is the electrostatic force acting on the unit charge placed at Fe of the heme. The 
requirement to reproduce electrostatic forces might be more stringent than to describe the elec- 
trostatic potential because of a shorter range and higher sensitivity to the local structure of the 
former. Since the statistical variances listed in Table [T] are not affected by the modeling of dissi- 
pative dynamics, the properties of the network itself need to be further fine-tuned. The inclusion 
of softening of the residue displacements by electrostatic water fluctuations (14) appears to be 
the first avenue for the model improvement. This addition will allow one to accommodate for 
the effect of the interface, in addition to the motions dictated by the shape and mass distribution 
which elastic models capture in the first place (|28l -l30h. Incorporating the water response will 
effectively produce a more heterogeneous distribution of the network force constants. 



CONCLUSIONS 

Network models of biomolecules are coarse-grained representations designed to describe their 
large-scale collective conformational motions (l27b . They capture the basic topology and pack- 



ing of residues in a biopolymer and robustly reproduce conformational global motions driven 
by the distribution of mass and shape (|28l-l30l). Electrostatic properties is yet another area where 



coarse-graining of biomolecules might be efficient. The long range of Coulomb interactions ef- 
fectively averages out details of the local structure suggesting that one can potentially describe 
the statistics and long-time dynamics of electrostatic fluctuations by global motions of charges 
assigned to the elastic network. The demonstration that this approach is in principle consistent 
with all-atom MD simulations is the main result of this study. Two components are critical to 
our approach: force-field atomic charges distributed at the residues of the network and two- 
exponential, overdamped dynamics assigned to each normal mode diagonalizing the network 
Hessian. 

Harmonic mechanical motions of the network of beads clearly do not incorporate dissipative 
dynamics of biomolecules in solution ([34, 45, 4^ and much is still needed to be done to achieve 



a physically consistent description of the protein dynamics. The elastic network employed here 
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assigns weaker springs between all non-covalent neighbors within a cutoff radius and stronger 
springs between covalent neighbors. Correspondingly, two global friction coefficients are as- 
signed to each normal mode diagonalizing the network Hessian. This phenomenological model 
qualitatively captures the two-peak loss spectrum of residue displacements and qualitatively 
similar loss spectra of the electrostatic potential and electric field. However, the characteristic 
adjustable friction coefficients used for the electrostatic fluctuations are not transferrable to the 
network displacements. A new set of parameters is needed when each property is considered 
separately. 
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